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ABSTRACT 


Moments  of  inertia  were  experimentally  determined  and  the  longitudinal  and 
lateral/directional  static  and  dynamic  stability  and  control  derivatives  were  estimated  for  a 
fixed  wing  Unmanned  Air  Vehicle  (UAV)  High  fidelity,  non-linear  equations  of  motion 
were  derived  and  tailored  for  use  on  the  specific  aircraft.  Computer  modeling  of  these 
resulting  equations  was  employed  both  in  Matlab/Simulink  and  in  Matrixx/Systembuild. 
The  resulting  computer  model  was  linearized  at  a  specific  flight  condition,  and  the 
dynamics  of  the  aircraft  were  predicted.  Several  flight  tests  were  conducted  at  a  nearby 
airfield  and  the  behavior  of  the  aircraft  was  compared  to  that  of  the  computer  model  The 
longitudinal  dynamics  as  depicted  by  the  short  period  mode  were  found  to  be  almost 
identical  with  those  predicted  by  the  nonlinear  computer  model.  The  phugoid  mode  was 
also  observed  and  found  to  be  in  close  agreement.  In  the  lateral/directional  dynamics, 
flight  test  was  employed  to  improve  the  model  and  the  parameters  were  modified  to  obtain 
a  better  match.  Ultimately  a  reasonably  accurate  nonlinear  model  was  achieved  as  required 
for  purposes  of  control  and  navigation  system  design 
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I.  INTRODUCTION 


A.  MISSION  NEED 

Unmanned  Aerial  Vehicles  (UAVs)  are  becoming  an  increasingly  important  part  of 
the  armed  forces  operations  During  various  military  missions  the  UAVs  have  the 
capability  to  collect  intelligence  and  target  for  gunfire  support,  gather  battle  damage 
assessment,  as  well  as  perform  many  other  tasks.  The  real  benefit  in  using  unmanned 
aircraft  lies  in  the  fact  that  no  lives  are  jeopardized  when  a  UAV  crashes  or  is  lost  to 
enemy  fire. 

Past  success  with  UAVs  in  the  Persian  Gulf  War  has  indicated  their  usefulness  in 
obtaining  timely  intelligence  during  military  conflicts  and  it  was  during  Desert  Storm  that 
the  first-ever  combat  tests  of  Unmanned  Aerial  Vehicles  that  were  obtained  by  U.S.  forces 
[Ref  1].  An  example  of  an  effective  UAV,  the  Pioneer  remotely  piloted  vehicle  (RPV), 
has  served  in  numerous  fleet  and  ground  operations  since  1987  [Ref.  2]. 

Another  advantage  of  UAVs  over  manned  aircraft  is  their  cost  which  is  only  a 
small  fraction  of  the  cost  of  a  manned  airplane.  Thus  they  have  become  an  integral  part  of 
modern  warfare  since  many  missions  such  as  electronic  deception,  visual  identification, 
laser  designation  of  targets  and  bomb  damage  assessment  deep  in  the  enemy  territory  can 
be  performed  without  endangering  any  lives  or  risking  any  expensive  aircraft.  They  can  be 
launched  from  practically  any  type  of  platform  making  them  ideal  for  use  in  the  Navy. 
They  are  also  usually  very  hard  to  be  detected  with  radar  or  infrared  systems  due  to  their 
small  size,  composite  materials  and  low  noise  and  speed.  UAVs  are  also  not  limited  by  the 
pilot  "g"  tolerance  or  fatigue. 


B.  EVOLUTION  OF  AN  AIRCRAFT 

For  any  type  of  aircraft  the  first  step  in  obtaining  an  accurate  mathematical  model 
is  to  determine  stability  and  control  derivatives.  These  derivatives  will  impact  the  flying 
characteristics  and  will  be  used  to  size  control  surfaces,  design  flight  control  systems  and 
program  devices  such  as  simulators. 

Three  approaches  can  be  taken  toward  completing  this  goal.  The  first  and  easiest 
method  requires  the  knowledge  of  the  geometry  and  inertial  properties  of  the  aircraft  and 
employs  simple  calculations  to  obtain  the  derivatives  within  reasonable  accuracy. 

The  second  method  involves  the  use  of  wind  tunnels.  However,  the  results  will 
have  to  be  refined  after  several  scale,  interference  and  dynamic  effects  are  taken  into 
account.  This  method  is  much  more  complicated  than  the  previous  one,  but  usually  more 
accurate. 

The  final  approach  which  is  the  most  time  consuming  and  costly  but  with  the  most 
promise  and  precision  is  flight  testing  Dynamic  flight  test  data  is  used  through  techniques 
such  as  parameter  estimation  to  accurately  estimate  the  stability  and  control  derivatives. 
Of  course  limitations  in  availability  of  data  and  noisy  measurements  can  cause  serious 
problems  in  the  successful  resolution  of  all  the  derivatives  of  interest. 

C.  REQUIREMENT  FOR  MODELING 

For  the  aeronautical  controls  engineer  the  goal  is  to  develop  a  dynamic  aircraft 
model  and  verify  its  accuracy.  This  model  will  then  be  used  to  design  a  control  system  for 
the  UAV. 

The  first  step  in  this  process  is  to  develop  the  high  fidelity  nonlinear  model  of  the 

aircraft  based  upon  the  predicted  flight  dynamics  involving  the  stability  and  control 

derivatives.  This  is  done  as  follows: 

•  The  equations  of  motion  are  developed  and  the  sum  of  all  forces  and  moments 
involved  are  obtained  allowing  for  an  easy  conversion  into  a  block  diagram 
representation. 


•  The  nonlinear  model  has  to  be  verified  through  flight  testing  and  if  necessary 
some  computer  aided  parameter  estimation  technique  should  be  incorporated  to 
obtain  accurate  results. 

Next  the  feedback  controller  is  designed  through  an  iterative  process  employing 

various  methods  with  the  design  requirements  of  response  time,  overshoot  and  command 

and  control  loop  bandwidths  being  the  adjusting  knobs  to  create  the  desired  controller. 

Flight  test  is  done  in  three  phases: 

•  On  a  laboratory  stand  allowing  only  a  restricted  degree  of  freedom  in  movement 
to  avoid  any  damage  to  the  aircraft. 

•  During  a  flight  with  an  experienced  pilot  ready  to  take  over  with  manual  control 
in  the  case  of  any  problem. 

•  Autonomous  flight  with  the  aircraft  in  full  operation. 

D.  STATEMENT  OF  OBJECTIVE 

Over  the  last  several  years  the  Avionics  Lab  of  the  Naval  Postgraduate  School  has 
undertaken  the  complicated  task  of  developing  and  implementing  control  systems  for 
UAVs.  Two  models  have  undergone  progress  up  to  the  point  of  conducting  flight  testing. 
The  first  one  was  the  Bluebird  which  was  a  conventional  aircraft  and  the  second  one  was 
the  Archytas  vehicle  with  a  ducted  fan  propulsion  system.  The  latest  vehicle  to  be  used  is 
the  Frog  unmanned  aerial  vehicle,  a  small  conventional  aircraft  acquired  as  a  test  bed  for 
designing  and  testing  guidance,  navigational  and  control  systems.  The  objective  of  this 
work  is  to  obtain  realistic  modeling  through  aerodynamic  parameters  of  the  aircraft  and 
through  the  acquisition  of  sufficient  flight  test  data  to  verify  the  nonlinear  model  which 
could  be  used  for  purposes  of  analysis  and  design  guidance,  navigation  and  control 
systems  Ultimately  these  guidance,  navigation  and  control  systems  along  with  a  Global 
Positioning  System  aided  autoland  capability  will  allow  for  fully  autonomous  operations. 


II.  AIRCRAFT  EQUATIONS  OF  MOTION  DEVELOPMENT 


A.    NOTATION 

First  it  is  necessary  to  present  the  notation  used  throughout  this  report  which  is 

consistent  with  the  previous  developments  [Refs.  6-8]  Consider  Fig.  2  1 

{A[  is  the  coordinate  system  with  basis  vectors  xa.va^a 

aPq  is  the  position  of  point  Q  expressed  in  J  A} 

AV(j  the  velocity  of  point  Q  measured  and  expressed  in  {A}. 

B(AV<j)  the  velocity  of  point  Q,  measured  in  {A}  and  expressed  in  {B} . 

BR  is  the  transformation  matrix  from  {B}  to  {A}. 

aQ.b  is  the  angular  velocity  of  the  (B|  coordinate  system  with  respect  to  {A} 
and  expressed  in  {A}. 

b(aQ.b)  is  the  angular  velocity  of  {B}  with  respect  to  {A}  but  now  expressed  in 


N^ 


{A} 


Ap 


,<•' 


BO 


Figure  2.1  Relative  Position  of  Coordinate  Systems 
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B.  COORDINATE  SYSTEMS 

The  following  coordinate  systems  and  assumptions  will  be  used  for  the  derivation 

of  the  equations  of  motion: 

{U}  is  the  tangent  plane  coordinate  system  attached  to  the  Earth. 

{B}  is  the  body  fixed  coordinate  system. 

{W}  is  the  wind  axis  coordinate  system. 

All  sensors  are  located  at  the  eg.  (This  can  be  relaxed;  however,  it  is  used  for 
simplification). 

The  mass  of  the  aircraft  remains  constant. 

For  the  vector  u  its  derivative  with  respect  to  {B}  is  jt(u)  and  with  respect  to. 
(U}is(«) 


The  {B}  coordinate  system  is  right-handed  with  X  aligned  with  the  thrust  axis.  All  rates  p, 

q,  r  are  defined  positive  when  clockwise  looking  in  the  positive  axis  direction.  The  {W} 

coordinate  system  has  X  aligned  with  the  wind  incident  on  the  aircraft.  Thus  a  is  the  angle 
formed  by  the  body  x-y  plane  and  the  positive  Xw  axis.  The  angle  |3  is  formed  by  the  body 
x-z  plane  and  the  positive  Yw  axis. 

For  simplification  the  following  definitions  are  introduced. 

ug  is  the  velocity  of  point  Q,  measured  and  expressed  in  {U} . 

v>bo  is  the  velocity  of  origin  of  {B},  measured  and  expressed  in  {U}. 

v)b  is  the  acceleration  of  {B}  with  respect  to  {U},  measured  and  expressed  in 
{U}. 

bvq  is  the  velocity  of  point  Q,  measured  in  {U}  and  expressed  in  {B} 

Ob  the  angular  velocity  of  {B},  measured  and  expressed  in  {U} . 

bg>b  is  the  angular  velocity  of  {B},  measured  in  {U}  and  expressed  in  {B} 

0  represents  the  appropriate  size  matrix  with  all  elements  zero. 

/„  is  the  identity  matrix  of  dimension  n. 


C.  EULER  ANGLES 

Using  the  three  Euler  angles  0, 0  and  *F  corresponding  to  roll,  pitch  and  yaw  the 
rotation  between  any  two   coordinate  systems  can  be  defined    For  our  case  of  a 

conventional  aircraft  the  3-2-1  order  of  rotation  is  used.  The  singularity  point  for  this 
order  arises  at  0  =  90°.  Thus  the  transformation  from  inertial  coordinates  {U}  to  body 

coordinates  {B}  can  be  carried  out  as  shown  in  Fig. 2. 2.  Therefore  the  direction  cosine 
matrix  is  expressed  in  terms  of  the  Euler  angles  as: 


cos  ¥  cos  0  sin  SP  cos  0  -sin0 

cos  *F  sin  0  sin  O  -  sin  ¥  cos  <t>  sin0sinOsinvF  +  coslFcos<I>  cos  0  sin  O 
cosxFsin0cosO  +  sinvFsinC>  sin 0 cos O sin *¥  -  cos *¥ sin O  cos0cosO 


(2.1) 


Next  the  kinematic  differential  equations  for  the  change  in  Euler  angles  are  given: 


P 

q 

r 


1       0  -sin© 

0   cos  O    cos  0  sin  O 
0  -sinO  cos©  cos  <5 


0 


(2.2) 


The  matrix  on  the  right  hand  side  is  invertible  for  all  0  *  nil  and  so  the  Euler  angle  rates 
can  be  obtained: 


0) 
0 


1  sinOtan©  cosOtan© 
0       cosO  -sinO 

0  sin O sec©  cosOsec© 


P 

q 

r 


(2.3) 


Therefore  the  time  history  of  the  Euler  Angles  is  obtained  integrating  equation  (2.3). 


Fig.2.2  Z-Y-X  Euler  Angles  Rotation  Sequence 


D.  DERIVATION  OF  EQUATIONS  OF  MOTION 

For  the  general  body  with  six  degrees  of  freedom  the  equations  of  motion  can  be 
derived  in  two  parts  In  the  first  part  is  the  determination  of  the  equations  of  motion  for  a 
rigid  body  and  only  the  linear  and  angular  momenta  is  considered  In  the  second  part  all 
the  forces  are  examined  and  more  specifically  aerodynamic,  gravitational,  and  thrust  forces 
are  taken  into  account.  It  is  in  the  modeling  of  the  aerodynamic  and  propulsive  forces 
where  the  stability  and  control  derivatives  come  into  picture 


1.  Linear  Equations 

These  are  derived  based  upon  the  Newton's  Law,  F=ma.  However  since  the 
measurements  are  resolved  in  the  body  {B}  coordinate  system  the  equations  are  also 
resolved  in  the  {B}  coordinate  system  Therefore,  using  Coriolis  theorem  we  obtain: 


BF=tnjt  \)B+mBG>BXBUB 


(2.4) 


2.  Angular  Equations 

To  derive  the  angular  momentum  equations  Euler's  Law  of  preservation  of 
momentum  is  employed.  Again  the  equations  are  written  in  the  {B}  body  coordinate 
system  for  the  position  of  the  eg.  of  the  vehicle.  Thus  we  get  [Ref  8]: 


3NBo  =  h  B®b+b  (Hb  x  (h  b&b) 


(2.5) 


3.  State  Equations 

Having  developed  the  equations  of  motion,  both  linear  and  angular,  the  next  step  is 
to  assemble  them  in  a  state  space  form  Thus  after  rearranging  and  normalizing  we  get: 


i 
(Ob 


-b  gsb*b  ubo 


-1  B 


+       ^ 
'  m 

1  B) 


rBl  bg)b  x  (iB  BaB)  +    rBl  bn 


(2.6) 


4.  External  Forces  And  Moments 

It  is  necessary  to  examine  now  the  external  forces  and  moments  acting  on  the 
rigid  body  and  distinguish  between  aerodynamic,  propulsive,  and  gravitational  forces: 


BF 

BN 


Bp-         ,B  p         ,B  f 

1  grav~     -»   prop'      1   a 

BN      +B  N 

**  prop^     •<' aero 


(2.7) 


a.  Gravitational  Forces 

These  act  on  the  body  and  have  to  be  rotated  by  the  appropriate  rotation 
matrix  from  the  tangent  to  the  body  coordinate  system.  This  matrix  is  defined  using  the 
Euler  angles: 


Br-         -B   D 
r  grav  —U-K 


0 

0 

mg 


(2.8) 


b.  Propulsive  Forces  And  Moments 

These  are  computed  directly  in  the  {B}  coordinate  system  and  are: 


Bf 
1  prop  — 


Tx 

Ty 


(2.9) 


and: 


3N 


prop 


Ti 

Tm 


(2.10) 


where  these  values  depend  on  the  aircraft  configuration  and  are  known. 

c  Aerodynamic  Forces  And  Moments 

The  derivation  of  these  forces  and  moments  depends  upon  the  nominal 

operating  point  about  which  they  have  to  be  found  using  a  Taylor  series  expansion.  The 

partial  derivatives  of  the  forces  and  moments  with  respect  to  the  aerodynamic  variables 
u/U,a,  fi,p,q,r  are  introduced  to  get  the  Taylor  series  expansion  [Ref  15]: 


Fa 


bFx>x'  +  hFx'x'  +  6FA  A  +  Fq 


(2.11) 
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Similarly  for  the  moments: 


Naero  =  §Nx/x'  +  hNx>x'  +  bNAA  +Nq 


(2.12) 


In  the  above  formulae  the  x'  is  the  state  vector: 


u 
P 

a 

pb_ 

2U 
qc_ 

2U 
rb_ 

2U 


(2.13) 


and  also: 


x'  = 


P 
d 


(2.14) 


Finally  for  the  control  vector: 


A  = 


8e 
6a 


(2.15) 


Combining  all  the  terms  together  and  introducing  the  matrix  of  non-dimensional  stability 
derivatives  differentiated  with  respect  to  the  terms  x',x',  A  we  get: 


§c_ 
dx' 


Cl»  Cl$  Clo.  ClP  CLq  Clt 

Cyv  Cyp  CVa  Cyp  Cyq  Cyr 

(--Do  Cz)p  Cz)a  (-Dp  (-Dq  (-Dr 

C/u  C/p  C/a  Clp  Clq  Clr 


Cr 


Cmf$     C 
C„p      C 


*--  rim    {--  oto    L  a 

^- np     t~-n<7     ^r 


(2.16) 
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dC 


The  term  —  is  similar  to  the  previous  term,  however,  only  the  terms  with  respect  to  a,  P 


are  computed  leaving  a  6x2  matrix.  Also  the  matrix  of  the  derivatives  with  respect  to  the 
control  inputs  elevator,  rudder  and  ailerons  is. 


ClZe     Cn>r     ClSq 

Cygg     Cyhr    Cyba 

dC  _ 

Code    Cobr    Co5a 

dA 

Clde      Cl8r      Cl8a 

L-  m8e    »--  mdr    ^  mba 

^  nhe     t>  rfor     ^  n&a 

(2.17) 


Finally  the  matrix  of  the  constant  coefficients  is: 


Cfo  — 


Cdo 
Cyo 
Clo 

Clo 

L  mo 


(2.18) 


The  above  values  refer  to  the  values  of  the  coefficients  at  the  trimming  point  and  not  to 
the  values  at  a  =  0°  [Ref.  10].  Also  the  stability  and  control  derivatives  are  found  in  the 
wind  axis  coordinate  system  and  it  is  necessary  to  transform  them  from  {W}  to  {B}.  The 

rotation  matrix  is  given  below: 


WR- 


cosa  cosP  -cos  a  sinP  -sin  a 

sin  P  cos  p  0 

sin  a  cosP    -sin  a  sinP    cos  a 


(2.19) 


The  aerodynamic  forces  and  moments  are  premultiplied  by  this  matrix.  In  addition  in  order 
to  be  consistent  with  the  way  lift  and  drag  are  defined  as  positive  we  have: 
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f  aero  — 


-D 

Y 
-L 


and    Nnero  = 


I 

m 
n 


(2.20) 


Finally  the  most  commonly  used  vector: 


x  = 


u 

V 

w 

p 
q 

r 


(2.21) 


is  introduced  and  the  complete  final  expressions  for  the  aerodynamic  forces  and  moments 
are  given: 


Bf 
r  a 


qS 


wR   o 
o  U 


{CF0+dC/dx' M'x  +  dC/dx' M'x  +  dC/dAA}  (2.22) 


where  M  and  M  map  from  x  to  x'  and  from  x  to  x'  respectively.  The  above  expression 
can  now  be  substituted  into  the  general  equation  2.6. 


5.  Complete  Equations  Of  Motion 

All  equations  presented  above  2.8,  2.9,  2.10,  2.22,  have  to  be  substituted  in  the 
general  equation  2.6  to  get  the  state  space  form.  After  introducing  the  terms: 


wT-- 


wR  o 
o  U 


and  Mr  = 


m    0 

0  BIB 


we  get  the  complete  set  of  equations  of  motion  which  in  the  state  space  form  are  the 
following  [Ref.  8]: 
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B 
©5 


X"1  {[ 


co5x 
0 


B  T-x(Bt 


0 

(BtB 


/iTfflax^/J©* +/>*)) 


+ 


Mjl*TqS%M'] 


VBO 

b 
(Ob 


+  MJl{ 


BF, 


grav 

0 


+ 


BF, 


prop 


'N. 


prop 


ecf 


$T+wTqS(CFo  +  ^-A)}] 


(2.23) 


Pbo  =r  R  bVbo 


A  =  S(A)bG)b 


X  =  h-M-jlBwTq~S^rM'. 


(2.24) 
(2.25) 
(2.26) 


where  P  is  the  position  vector  of  the  aircraft,  and  5(A)  is  the  matrix  of  kinematic 
differential  equations  based  upon  the  Euler  angles. 
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III.  THE  AIRCRAFT 


A.  GENERAL  DESCRIPTION 

The  Frog  unmanned  aerial  vehicle,  shown  in  Figure  3.1,  is  a  high-wing 
tricycle-gear  radio-controlled  airplane.  It  is  constructed  of  wood,  foam,  composites  and 
metal.  It  is  powered  by  a  two-stroke,  air-cooled  engine  with  a  shaft  horsepower  of  5.6  hp. 
It  is  controlled  by  a  Futaba  flight  control  transmitter  operating  on  a  frequency  of  72  MHz 
To  enhance  reliability,  a  factory  variant  of  the  control  system  is  available  that  provides 
direct  control  from  a  transmitter  that  transmits  on  407.275  MHz  and  this  configuration 
could  be  used  in  the  event  that  the  primary  relay  should  fail  when  the  aircraft  is  over  the 
horizon.  Table  3.1  describes  the  physical  characteristics  of  the  Frog. 


Length 

8.125  ft 

Height 

1.75  ft 

Wing  Airfoil  (est.) 

NACA2415 

Horizontal  Stab.  Airfoil  (est.) 

NAC  A  0010-34 

SwingWref) 

17.5  ft2 

s. 

3.175  ft2 

Sv 

0.9818  ft2 

c 

1.66  ft 

c, 

0.968  ft 

b 

10.58  ft 

b, 

3.313  ft 
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bv 

1.25  ft 

AR 

6.37 

AR, 

3.46 

ARV 

1.59 

vH 

0.485 

Vv 

0.0235 

V 

0.0022 

1, 

4.44  ft 

lv 

4.44  ft 

Table  3.1.  Specifications 

B.  STABILITY  DERIVATIVES 

The  flight  condition  for  which  the  aircraft  model  was  obtained  is  described  in  Table 
3.2.  Based  upon  these  values  the  initial  estimates  of  the  stability  derivatives  were  made 
using  the  physical  characteristics  of  the  aircraft  such  as  airfoil  data,  geometric 
measurements,  relative  positions  of  aircraft  components,  mass  and  weight  [Refs.  11-14]. 
Theoretical  or  quasi-theoretical  methods  taken  from  the  literature  were  used  for 
calculating  these  constants  thus  forming  the  basis  for  the  first  parts  of  longitudinal  and 
lateral  stability  determinations  [Ref  12].  These  methods  are  regarded  as  the  most  suitable 
for  light  aircraft  with  a  typical  configuration. 


Weight 

67.73  lbs. 
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^xx 

12.52  slug*ft2 

*YY 

8.43  slug*ft2 

*ZZ 

18.55  slug*ft2 

Velocity 

60  mph/88  f7sec 

Altitude 

800  ft  MSL 

Air  Density 

0.002327  slugs/ft3 

Center  Of  Gravity 

34.5%  of  mac. 

c 

0.4295 

A  0  Atrim 

5.25° 

Elevator^ 

5.14° 

Table  3.2  Flight  Condition/ Aircraft  Configuration 

Matlab  programs  [Ref  21]  written  for  the  physical  and  derivative  calculations  are 
presented  in  Appendix  A  Next  the  nondimensional  stability  and  control  derivatives  are 
shown  in  Table  3.3  and  dimensional  stability  and  control  derivatives  estimated  are  shown 
in  Table  3.4. 
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cL 

0.2866 

cD 

0.0614 

Cu 

4.3034 

Coo 

0.0499 

c 

1.3877 

c 

-3.7115 

Ladot 

Madot 

^Ma 

-0.5565 

CD. 

0.23 

Cl, 

3.35 

^Mq 

-8.8818 

^Lde 

0.3914 

c 

0.0676 

Dde 

C     , 

-1.0469 

c* 

-0.31 

mde 

C* 

0.0575 

c 

-0.0509 

Cyp 

0.0000 

c* 

-0.0537 

c 

^1P 

-0.3702 

c. 

0.1151 

cm 

-0.0669 

c„ 

0.1119 

Cyda 

0.0000 

^nda 

-0.0272 

Cy. 

0.1810 

^ydr 

0.0926 

^ndr 

-0.0388 

cldr 

0.0036 

cDq 

0.0000 

c* 

0 

Table  3.3  Nondimensional  derivatives 
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x„ 

-0.1045  /sec 

x4 

14.9485  ft/sec 

xde 

-5.0602  ft/sec2 

zu 

-0.7312  /sec 

z. 

-326.9292  ft/sec2 

z,t 

0.9804  ft/sec 

»dot 

z, 

2.3667  ft/sec 

Zde 

-29.3185  ft/sec2 

H, 

0.0000  /ft*sec 

M, 

-17.2801  /sec2 

M„, 

-1.0869 /sec 

M, 

-2.6010 /sec 

»dot 

Mde 

-32.5049 /sec2 

YB 

-23.2196  ft/sec2 

YP 

0.0000  ft/sec 

Yr 

0.5181  ft/sec 

Ya. 

0.0000  ft/sec2 

Y* 

-6.9323  ft/sec2 

Lb 

-6.7831  /sec2 

Lp 

-2.9653  /sec 

Lr 

0.8964 /sec 

L,u 

24. 1120 /sec2 

L* 

0.4849  /sec2 

NB 

5. 1744 /sec2 

NP 

-0.2903  /sec 

Nr 

-0.3619 /sec 

N* 

-2.4467  /sec2 

N* 

-3.4927 /sec2 

Table  3.4  Dimensional  derivatives 
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Fig.  3.1  Frog  Drawing 
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C.  MOMENTS  OF  INERTIA 

When  developing  a  mathematical  model  of  the  UAV  it  is  critical  to  have  accurate 
moments  of  inertia.  Direct  calculation  of  a  model's  moments  of  inertia  by  consideration  of 
the  contributions  made  by  individual  parts  is  rather  impractical  and  inaccurate  because  the 
model's  parts  are  too  small  and  light  and  the  distances  too  short  to  yield  anywhere  near 
accurate  values  for  the  moments  of  inertia  [Ref  5].  That  is  why  testing  is  employed  to 
determine  the  moments  of  inertia  more  precisely  in  practice.  Any  changes  that  might  occur 
in  a  model's  moment  of  inertia  due  to  addition  or  substraction  of  equipment  or  structure 
could  be  calculated  directly  thereafter 

Moment  of  inertia  is  a  measure  of  a  rotating  body's  resistance  to  acceleration  and 
can  be  computed  by  taking  the  product  of  the  mass  of  each  part  of  the  aircraft  and  the 
square  of  the  distance  from  the  aircraft's  center  of  gravity.  There  are  two  approaches  to 
obtaining  the  moments  experimentally,  each  employing  the  principle  of  the  compound 
pendulum  and  each  giving  the  same  answer.  One  method  involves  hanging  the  model  from 
a  single  point  in  the  ceiling  on  two  wires  or  cords,  one  well  forward  of  the  center  of 
gravity  and  the  other  well  aft.  Thus  a  compound  pendulum  is  obtained.  Another  technique 
is  to  make  the  center  of  gravity  of  the  model  itself  the  pivot  point  and  complete  the 
compound  pendulum  by  hanging  a  bob  weight  below  it  on  a  pair  of  fairly  stiff  wires  or 
lightweight  (wooden)  struts.  Needless  to  say  that  the  friction  at  the  pivot  point  should  be 
held  to  a  minimum.  Next  by  giving  a  small,  gentle  push  in  the  appropriate  direction  and 
timing  its  oscillations,  the  oscillatory  period  (T)  is  determined  which  is  simply  the  total 
number  of  seconds  divided  by  the  total  number  of  cycles  (one  cycle  is  one  complete  swing, 
to  and  fro).  It  should  be  mentioned  that  the  greater  the  number  of  cycles  (at  least  20-30) 
and  the  longer  the  suspension  the  greater  the  timing  accuracy,  which  is  vital. 

Using  the  period  exhibited  and  the  principles  described  above  the  moments  of 
inertia  were  calculated  [Refs.  5,  22].  In  order  to  calculate  the  moments  about  all  three 
axes,  the  model   was  hung  and  swung  three  different  ways,  each  time  about  the  axis  of 


21 


interest.    It   was   hung   by   chain   and   swung   as   pictured   in   Figures  3.2,    3.3,   3.4. 
Specifications  for  the  geometry  of  each  test  can  be  found  in  Appendix  B 

Reference  23  provides  equation  3.1  for  calculating  the  moment  of  inertia  for  a 
swinging  model: 

/  =  ^~ -P^—zr — —  (31) 

where  W  is  the  weight,  Z  is  the  distance  from  the  pivot  point  to  the  center  of  gravity,  p  is 
the  period,  and  g  the  gravitational  constant.  M,  S  are  subscripts  referring  either  to  the 

model  or  the  support. 

It  was  determined  that  swinging  just  the  support  (chains),  in  the  configuration  it 
would  be  in  when  supporting  the  model  would  not  be  possible,  since  the  chains  would  not 
maintain  their  positions  without  the  model  in  place.  Equation  (3.1)  was  modified  in  order 
to  treat  the  chains  as  long  slender  rods  and  to  calculate  their  moments  of  inertia  as  such 
[Ref.  22].  The  new  form  of  equation  is: 


I=-^—Pu+s--g--Z^r-  (3.2) 

where  Ls  is  the  length  of  the  chain  and  the  summation  is  taken  over  all  chains  (four  in  this 
case).  In  particular,  two  long  and  two  short  chains,  all  with  a  weight  per  unit  length  o. 

After  all  the  appropriate  substitutions  were  made: 


J  _   [WM±2<£>(LsHORT+LiX)NG)]Z.MJrS  ^2  WM%M  2a>  /  j  3  ,    J  3  \  /"J   -2\ 

y-  ^2  PM+S  g  3JV^  SHORT   +L'LONG)  \A-*) 

Having  the  equation  in  this  form  all  variables  are  fixed  except  Zm+s,  Zm,  Pm+s  Therefore 
for  each  configuration  these  three  variables  were  measured  and  all  necessary  calculations 

were  made.  Three  periods  were  timed  during  the  swing  tests  and  the  tests  were  repeated 

ten  times.  The  moments  of  inertia  calculated  are  shown  in  Table  3  2 
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Fig.  3.2  1^  Test 
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Fig.  3.3  I^Test 
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Fig.  3.4  lm  Test 
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IV.  COMPUTER  MODELING 


In  Chapter  II  the  full  nonlinear  equations  of  motion  were  developed  and  in  Chapter 
III  the  complete  set  of  stability  and  control  derivatives  were  obtained  Now  the  next  step 
is  to  develop  the  model  of  the  aircraft  on  a  computer.  This  model  should  be  in  a  generic 
format  and  should  accept  values  for  derivatives  from  a  generic  input  file  Thus  by  simply 
changing  the  input  data  for  any  aircraft  the  model  could  be  used  for  a  different  aircraft. 

For  this  report  two  already  existing  models  was  employed,  the  first  one  developed 
by  DR.  Kuechenmeister  implementing  the  equations  of  motion  in  Matlab  and  Simulink 
[Ref.  6]  and  the  second  one  in  Systembuild  [Ref.  19].  These  models  have  been  validated 
by  inputting  the  appropriate  data  for  a  well  known  aircraft,  such  as  the  Cessna  172  and 
comparing  the  results  of  the  model  to  existing  data.  It  was  found  that  the  results  from  the 
numerical  linearization  are  quite  close  and  therefore  the  equations  used  in  the  model  are 
considered  to  be  correct  [Ref.  6]. 

A.  BASIC  NONLINEAR  MODEL 

1.  Basic  Simulink  Model 

The  basic  nonlinear  model  was  constructed  by  implementing  the  state  derivative 
equations  in  SIMULINK.  This  model  is  shown  in  Figures  4. 1  and  4.2  below.  In  this  model 
the  Matlab  function  blocks  abmat.m  and  frogdata.m  are  used  to  input  the  equations  of 
motion  and  the  stability  and  control  derivatives,  respectively.  Notice  that  the  force  due  to 
thrust  since  no  previous  measurements  existed,  was  assumed  to  be  equal  to  drag.  It  was 
decided  to  follow  this  simplified  approach  at  this  stage  since  no  sufficient  information  was 
given  about  the  engine  in  terms  of  torque  or  propeller  parameters.  In  the  previous  works 
the  term  of  thrust  was  computed  from  an  assumed  propeller  efficiency  and  the  given 
horsepower  of  the  engine  [Refs.  6,  7]  However  it  has  been  found  [Ref.  4]  that  especially 
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the  value  of  horsepower  of  the  engine  given  from  the  manufacturer  is  far  from  accurate 
and  is  not  attained  in  flight.  A  more  detailed  description  of  the  procedure  for  determining 
the  shaft  power  from  ground  torque  tests  and  wind  tunnel  tests  for  propeller 
characteristics  is  given  in  [Ref.  4].  Therefore  the  propulsive  forces  were: 


Fprop  - 


D 
0 
0 


ST 


(4.1) 


and  for  the  propulsive  moments  using  the  position  of  the  propeller  from  the  center  of 
gravity  we  obtained  that: 


NpROP^PenfD-ST 


(4.2) 


where  BPeng  is  the  position  of  the  engine  in  {B}  body  coordinate  system.  In  Appendix  C 
all  the  input  values  are  found  in  file  Frogdat  m  and  the  file  Abmat.m  for  the  equations  of 

motion. 


Fig.  4.1  Equations  of  Motion  Implementation 


28 


Next  the  Simulink  diagram  for  the  block  of  nonlinear  equations  of  motion  of  Figure  4.1  is 
presented  in  Figure  4.2.  In  this  model  the  airspeed  and  the  flight  path  angle  have  been 
selected  as  outputs  since  they  will  be  used  later  in  the  trimming  and  linearization  of  the 
model. 


Q- 


& 
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Fig.  4.2  Nonlinear  Equations  of  Motion  Block 


2.  Basic  Systembuild  Model 

The  state  derivative  equations  were  also  implemented  on  Systembuild  in  a  similar 
manner.  Xmath/Systembuild  is  a  software  program  similar  to  the  Matlab/Simulink 
software  program  developed  by  Math  Works  Inc.  It  is  suitable  for  both  the  classical 
input/output  control  techniques  and  the  modern  state-space  representations.  Specifics  can 
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be  found  in  [Ref  19]  and  the  Xmath  and  Systembuild  Core  manuals.  With  Systembuild 
using  a  hierarchical  method  of  organization,  based  on  the  Superblock  concept  any  model 
can  be  drawn  and  tested.  Therefore  the  nonlinear  equations  of  motion  were  implemented 
as  shown  in  Figures  4.3.  All  the  Superblock  diagrams  are  presented  in  Appendix  E. 

In  the  superblock  of  the  Integrators  the  derivatives  of  the  state  vector  are 
integrated  to  obtain  linear,  angular  velocity,  the  Euler  angles  and  position  in  cartesian 
coordinates.  The  dynamics  are  implemented  in  the  superblock  "DynamicsEuler"  and  the 
data  for  the  specific  aircraft  is  input  to  the  superblock  of  "Aeroforcesandmoments" 
inside  the  superblock  of  "DynamicsEuler".  Thus  the  differential  equations  of  the  linear, 
angular  and  Euler  angles  equations  are  obtained  in  the  corresponding  superblocks  which 
produce  the  derivatives  of  the  state  vector.  Moreover  the  same  assumptions  concerning 
the  thrust  applied  by  the  engine  were  made  as  in  the  Matlab/Simulink  model 
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Fig.  4.3  Frog  Systembuild  diagram 
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B.  LINEARIZATION  OF  DEVELOPED  MODEL 

The  linearization  process  started  with  the  trimming  of  the  equations  of  motion  at 

the  nominal  flight  condition  of  VT  =  88/// sec  and  y  =  0.  Furthermore  it  was  also  selected 
that  for  a  cruise  condition  the  bank  angle  should  be  (j)  =  0  and  sideslip  (3  =  0  Therefore  the 

complete  set  of  the  nine  equations  of  motion  were  solved  for  the  rest  of  the  state  vector 
and  input  vector  which  were  unknown  in  trim  using  the  "trim"  command  in  Matlab: 

x0  =  [  88  0  0.1593  0  0  0  0  0.0018  0  ]' 

i/0  =  [  -0.0431   0  0  0.9805  ]' 

then  the  output  will  be  (airspeed/sideslip): 

yo  =  [88,or 

Next  the  linearization  was  obtained  in  Matlab  using  the  "linmod"  command  and  the 
matrices  A,  B,  C,  D  were  found  along  with  the  eigenvalues  of  the  A  matrix  and  are 
presented  in  Table  4.1. 


Eigenvalues 

Damping 

Freq.(radVsec) 

0.0474 

-1.0 

0.0474 

0 

-1.0 

0 

-0.0293+0.5597i 

0.0522 

0.5604 

-0.0293-0.5597i 

0.0522 

0.5604 

-0.1964+2.4972i 

0.0784 

2.5049 

-0.1964-2.4972i 

0.0784 

2.5049 

-3.2691 

1.0 

3.2691 

32 


-3. 7591+3. 5964i 

0.7226 

5.2024 

-3.7591-3.5964i 

0.7226 

5.2024 

Table  4.1  Eigenvalues  of  Linearized  model 

The  complete  numerical  results  of  the  linearization  process  are  presented  in  Appendix  D. 
It  is  easy  to  identify  the  various  modes  of  the  longitudinal  and  lateral/directional  dynamics 
of  the  aircraft  [Ref  24],  The  two  pairs  of  complex  eigenvalues: 

Xsp  =  -3.7591  ±3.5964/ 


and: 


XP  =  -0.0293  ±0.5597/ 

correspond  to  the  short  period  and  phugoid  modes  respectively.  It  may  be  noted  that  the 
short  period  is  very  highly  damped  and  with  a  natural  frequency  of  5.2024  radl sec.  It  is 
therefore  within  the  satisfactory  boundaries  of  the  short  period  "thumbprint",  although  the 

response  initially  could  be  a  bit  abrupt.  The  phugoid  mode  is  very  lightly  damped 
(C  =  0.0522)  as  expected  and  with  a  very  low  natural  frequency  of  0.5604  radl  sec . 

In  the  lateral/directional  dynamics  the  dutch  roll  mode  can  be  easily  identified  with 

the  following  pair  of  complex  eigenvalues: 

k  =  -0.1 964  ±2.4972/ 
and  has  a  damping  ratio  and  natural  frequency: 

Qd-r  =  0.0784,        On  =  2.5049  radl  sec 


The  dutch  roll  damping  is  characterized  as  low  to  moderate  and  therefore  the  response  is 
apparent  but  should  not  give  serious  difficulty  in  maneuvers.  The  value  of  the  natural 
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frequency  in  dutch  roll  is  moderate  to  high  and  in  some  cases  the  airplane  may  become 
oversensitive. 

We  also  observe  that  the  roll  response  has  a  stable  eigenvalue  of: 

X  =  -3.2691 

and  therefore  has  a  time  constant  of : 

T=  1/^  =  0.3059  sec 

for  time  to  half  amplitude 

Finally  the  only  unstable  mode  is  the  spiral  with  an  eigenvalue: 

^spiral  =  0.0474 

with  a  time  constant  of  21 .097  sec  which  is  considered  as  rather  low  Due  to  a  rather  large 
derivative  of  Np  the  spiral  mode  is  divergent  and  may  affect  the  flying  qualities  of  the 
aircraft  since  it  could  result  in  difficulty  in  lateral  trimming  for  wings  level  flight. 
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V.  FLIGHT  TESTING 


The  primary  goal  of  this  research  was  to  investigate  the  flying  qualities  of  the  UAV 
Frog.  The  nonlinear  model  which  has  been  developed  cannot  be  used  to  design  any  useful 
and  realistic  controller  unless  flight  testing  verified  its  accuracy.  It  is  therefore  desirable: 

•  To  correlate  the  analytic  parameter  estimates  with  flight  test  data. 

•  To  more  accurately  refine  the  parameter  estimates  for  purposes  of  control 
system  analysis  and  design. 

It  seems  wise,  therefore,  to  use  flight-test  data,  at  the  very  least,  as  a  verification  tool  of 

aircraft  stability  and  control  derivatives.  Ultimately  a  more  sophisticated  computer-aided 

parameter  estimation  routine  will  be  incorporated  that  will  lead  to  the  accurate  extraction 

of  the  stability  and  control  derivatives. 

Flight  tests  were  conducted  using  the  UAV  Frog  at  an  outlying  field  near  the 

Naval    Postgraduate    School.    These   tests   involved   transporting   the    Sun    Sparc    2 

workstation  Intrepid,  the  luggable  PC  AC  100,  the  communication  box,  the  RF  antenna, 

the  portable  generator,  and  the  airplane  to  the  field.  All  the  appropriate  connections  and 

necessary  steps  are  described  in  [Ref  19].  Before  any  flight  test  could  be  conducted,  all 

the  calibrations  of  actuators  and  sensors  were  ensured.  The  actual  flight  test  was  then 

commenced  and  data  collection  was  initiated.  The  data  saved  was  finally  converted  to  a 

format  suitable  to  be  analyzed  in  Xmath.  One  of  the  benefits  of  the  Matrix  software  was 

the  ability  to  collect  and  analyze  flight  test  data  during  the  flight  testing  process.  For  a 

complete  description  of  all  steps  see  [Ref.  19]. 

A.  LONGITUDINAL  DYNAMICS 

Due  to  limitations  in  time  and  resources  one  dynamic  mode  of  the  longitudinal 
dynamics  was  mostly  studied  The  short  period  mode  of  an  airplane  is  the  one  of  the  two 
modes  of  the  longitudinal  dynamics  and  concerns  the  pitching  motion  of  the  plane  about 
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the  center  of  gravity  with  little  or  no  airspeed  variation,  characterizing  the  ability  to  return 
to  the  trim  angle  of  attack  following  some  disturbance.  The  other  longitudinal  mode,  the 
long  period  mode  or  phugoid  is  a  gradual  interchange  between  potential  and  kinetic 
energy.  The  phugoid  is  in  general  slow  and  controllable  by  the  pilot  while  the  short  period 
is  not  and  therefore  subject  to  tight  specifications.  Since  the  short  period  is  the  mode 
affecting  most  of  the  longitudinal  flying  qualities,  it  was  singled  out  for  study. 

To  evaluate  the  aircraft  performance  the  instrumentation  measured  the  following 
critical  parameters: 

♦  angle  of  attack,  a 

♦  pitch  rate,  q 

♦  elevator  deflection,  5e 

Also  the  airspeed  could  only  be  obtained  through  the  GPS  data  recordings  since  no 
accurate  airspeed  measurement  unit  was  installed  when  the  tests  took  place.  It  has  also 
been  found  that  the  pendulums  used  to  measure  the  angles  give  accurate  results  only 
during  steady  state  conditions  and  could  not  be  relied  during  any  kind  of  maneuver, 
therefore  their  measurements  are  not  presented. 

The  measured  and  recorded  elevator  deflection  was  the  input  to  the  nonlinear 
model  of  the  aircraft  in  Matrix,^  and  the  outputs  consisting  of  angle  of  attack  and  pitch  rate 
were  plotted  versus  the  flight  test  results.  Five  individual  runs  were  conducted  during  the 
flight  test,  and  for  each  run  an  elevator  doublet  was  applied  [Ref.24].  To  excite  the  short 
period  mode  only,  the  goal  was  to  have  as  little  deviation  as  possible  in  pitch  attitude  from 
trim  at  the  end  of  the  doublet.  The  technique  used  was  as  follows: 

♦  trim 

♦  apply  5-10°  nose  down,  then  apply  5-10°  nose  up 

♦  release  and  record  the  aircraft's  response 

In  Figure  5.1  the  results  for  the  first  run  are  presented.  The  full  results  for  all  runs  can  be 
found  in  Appendix  F. 
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Fig.  5.1  Nonlinear  model  response  plotted  over  flight  test  data/  First  maneuver 
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It  is  apparent  that  there  is  a  close  agreement  between  the  flight  test  data  and  the 

response  of  the  nonlinear  model.  In  both  cases  the  short  period  damping  was  quite  high, 

resulting  in  an  absence  of  overshoot  while  the  initial  response  is  very  slow  since  the  natural 

frequency  has  a  moderate  value  [Ref.  24].  (In  some  cases  the  peaks  of  the  aircraft's 

response  in  both  angle  of  attack  and  pitch  rate  look  like  they  have  been  cut  off  and  this 

could  be  identified  as  mostly  a  sensor  problem )  As  another  aside  a  constant  bias  of 
OLbias  =  15°  in  angle  of  attack  was  added  in  all  cases  due  to  a  difference  in  trim  values. 
Nevertheless  the  close  agreement  verifies  the  values  of  the  longitudinal  parameters 

calculated  previously  analytically.  It  is  therefore,  a  valid  nonlinear  longitudinal  model  for 

purposes  of  control  system  analysis  and  design  [Ref.  25]. 

Although  the  short  period  mode  was  mostly  pursued,  for  comparison  reasons 
mainly,  the  response  of  the  aircraft  in  airspeed  and  altitude  was  plotted  versus  the  one 
obtained  from  the  nonlinear  model.  The  GPS  data  was  employed  since  no  accurate 
information  was  obtained  for  airspeed  and  altitude  from  the  Inertial  Measurement  Unit. 
From  the  Systembuild  model  it  was  easy  to  extract  the  airspeed  and  altitude.  Due  to  the 
much  lower  frequency  the  GPS  data  exhibits  a  step  behavior.  In  Figure  5.2  the  altitude  and 
airspeed  variations  are  plotted  for  both  the  flight  test  and  the  simulation  for  the  first 
maneuver,  while  the  full  set  of  plots  are  presented  in  Appendix  F.  It  was  found  that  they 
are  in  a  relatively  close  agreement  and  thus  the  longitudinal  nonlinear  model  is  verified  for 
the  phugoid  mode  as  well. 

For  even  more  accurate  predictions  of  the  longitudinal  control  and  stability 
derivatives  a  computer  aided  parameter  estimation  program  should  be  incorporated  once  a 
full  state  vector  could  be  measured  in  flight.  This  will  provide  the  ability  to  estimate  the 
parameters  in  the  presence  of  state  noise;  however,  a  more  accurate  and  comprehensive 
data  acquisition  system  should  be  employed  for  this  purpose  [Refs.  17,  18] 
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Fig.  5.2  Phugoid  response  of  nonlinear  model  over  flight  test  data/  First  maneuver 
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B.  LATERAL  DIRECTIONAL  DYNAMICS 

1.  Steady  Heading  Sideslips 

To  study  the  behavior  of  the  aircraft  in  terms  of  the  lateral/directional  dynamics  as 
well  as  to  evaluate  the  fidelity  of  the  nonlinear  computer  model  a  series  of  steady  heading 
sideslip  maneuvers  were  conducted.  By  measuring  the  control  deflections  and  forces  that 
are  required  to  hold  the  aircraft  away  from  trim  the  equal  and  opposite  forces  that  tend  to 
stabilize  the  plane  can  be  found.  During  this  steady  maneuver  the  yaw  and  roll  rates  are 
both  zero  and  therefore  their  contributions  are  also  zero.  The  relations  governing  this 
maneuver  are: 

CrBP  +  Cy&rdr  +  CY5ada  =  -CL  ■  sin(4>)  (5.1) 

C„pp  +  C„8r5r  +  CBSfl  =  0  (5.2) 

C/„P  +  C/8r8r  +  C/8.=0  (5.3) 

During  the  flight  test  the  following  parameters  were  recorded: 

•  sideslip,  p 

•  airspeed 

•  rudder  deflection,  5r 

•  aileron  deflection,  80 

Using  the  measured  values  of  rudder,  aileron  deflections  and  bank  angle    for  various 

values  of  sideslip,  the  plots  of  deflection,  bank  angle  versus  sideslip  were  obtained  in 

Figure  5.3.  The  response  was  found  to  be  fairly  linear  as  expected  and  the  slopes  of  these 

lines  were  extracted.   The  same  slopes  were  also  computed  from  the  theoretically 

calculated  values  of  the  derivatives  using  the  above  equations  5.1-5.3  and  found  to  be  in 
quite  close  agreement  with  the  exception  of  the  slope  <J)/f3: 
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Ratio 

Flight  Test 

Theoretical 

&V(3 

1.286 

1.303 

Sa/(3 

0.265 

0.255 

w 

1.353 

0.441 

Table  5.1-  Sideslip  Results 

Based  upon  the  results  for  the  slopes  from  the  flight  testing  the  values  of  the  (3 
derivatives  were  calculated  from  equations  5.1-5.3.  This  yielded  improved  values  of  the 
Crp,  C/„,  C„9  derivatives:  -0.70017,  -0.05254  and  0.057085  respectively  and  in  Table 
5.2  the  old  and  the  improved  values  are  presented.  With  these  values  the  steady  heading 
sideslip  maneuvers  were  simulated  in  the  nonlinear  model  and  the  response  was  compared 
to  the  actual  aircraft  behavior  in  Figure  5.4  while  the  full  set  of  these  plotted  responses  can 
be  found  in  Appendix  F.  It  is  easily  seen  that  the  simulation  is  in  quite  good  agreement 
with  the  flight  test  data.  Gust  effects  are  not  modelled  and  may  provide  some  discrepancy 
in  the  analysis.  It  is  also  apparent  that  the  level  of  electronic  noise  is  quite  significant  in 
both  the  response  of  the  aircraft  and  the  control  inputs  and  to  avoid  unrealistic  inputs  to 
the  model  the  control  input  data  was  processed  by  deleting  any  spikes  caused  by  noise. 
However  this  is  a  steady  state  response  and  in  general  there  is  an  acceptable  match 
between  the  nonlinear  model  and  the  plane's  response. 


Derivative 

Old  value 

New  value 

CV(5 

-0.31 

-0.70017 

C/p 

-0.0509 

-0.05254 

t-np 

0.0575 

0.057085 

Table  5.2-  Sideslip  derivatives  comparison 
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Fig.  5.3  Rudder,  Aileron,  Bank  angle  vs.  Sideslip  Angle,  trimmed  sideslips 
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Fig.  5.4  Steady  Sideslip  Response  of  Computer  Model  over  Flight  Test  Data 
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2.  Dutch  Roll/Spiral  Response 

To  further  investigate  the  fidelity  of  the  computer  model  the  dutch  roll  mode  was 

excited  in  flight  testing  by  applying  a  rudder  doublet  and  after  the  plane  settled  for  some 

time  applying  an  aileron  doublet,  leading  to  the  spiral  response.  The  same  response  was 

obtained  from  the  model  and  the  results  in  terms  of  the  roll  rate  and  yaw  rate  are 

presented  in  Figure  5.5  and  in  Appendix  F.  It  was  found  that  adjustment  of  the  damping 
derivatives:  Cip,  Cir,  C^and  C„r  was  necessary  to  get  a  relatively  close  match  between 
the  responses.  Initially  the  computer  model  was  much  less  damped  than  the  actual  aircraft 

and  with  a  different  period.  By  adjusting  these  derivatives  through  an  iterative  process  the 

behavior  of  the  computer  model  was  improved,  however,  further  processing  is  required  to 

achieve  a  better  agreement.  The  response  in  terms  of  the  bank  angle  is  not  shown  since  the 

pendulums  used  in  the  inertial  navigation  unit  are  not  trustworthy.  Another  problem  was 
that  no  channel  was  available  during  this  maneuver  to  record  the  sideslip  P  since  this  is  not 
obtained  as  one  of  the  IMU  parameters.  After  adjusting,  the  improved  values  of  the 
damping  derivativesC/,,  C/r,  C^and  C„rare:  -0.3702,  0.4476,  -0.1074  and  -0.1338 
respectively.  It  was  found  that  it  was  necessary  to  decrease  them  (more  negative)  in  order 
to  make  the  nonlinear  model  more  damped  and  match  it  with  the  actual  plane's  response 


Derivatives 

Old  values 

New  values 

C; 

-0.3702 

-0.3702 

Ci, 

0.1119 

0.4476 

Cn, 

-0.0537 

-0.1074 

c„. 

-0.0669 

-0.1338 

Table  5.3  Damping  Derivatives  Comparison 
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Fig.  5.5  Angle  Rates  Response  in  Rudder/ Aileron  Doublet 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

A  determination  of  the  moments  of  inertia  through  the  compound  pendulum 
analysis  was  done  with  reasonable  results  Future  changes  to  the  configuration  of  the 
aircraft  could  be  accounted  for  by  adjusting  for  the  contribution  of  each  individual 
component  added  or  removed 

Initial  estimates  of  stability  and  control  derivatives  were  made  by  means  of 
conventional  aircraft-design-type  methods.  Taking  into  account  characteristics  such  as 
lift-curve  slopes,  geometry  and  weight  parameters  measured  the  derivatives  were 
estimated  in  Matlab  with  programs  which  allow  for  any  future  changes  or  for  considering 
different  flight  conditions. 

A  high  fidelity  nonlinear  model  of  the  Frog  was  implemented  in  Matrixx 
/Systembuild  using  the  superblock  diagram  structure.  This  structure  allows  for  changes  to 
be  easily  made  and  requires  little  or  no  programming  ability,  other  than  some  familiarity 
with  this  software.  One  significant  benefit  of  this  product  is  the  ability  to  collect  and 
analyze  data  even  during  the  flight  testing  process.  This  allowed  us  to  make  changes  and 
record  different  data  in  the  field  without  having  to  dismantle  the  equipment  The  real-time 
data  acquisition  allowed  to  convert  the  data  of  the  flight  test  to  a  form  suitable  for  analysis 
purposes.  Since  all  inputs  and  outputs  could  be  recorded  at  each  time  step  thus  the  data 
was  scrutinized  after  each  test. 

Longitudinal  flight  tests  supported  with  a  great  success  the  fidelity  of  the 
parameters  estimated  both  in  short  period  and  in  phugoid  modes.  The  results  of  the 
simulation  were  almost  identical  to  the  flight  testing  data  and  were  both  dictated  by  the 
strong  damping  and  the  limited  response.  In  the  lateral-directional  tests  the  comparison 
between  the  simulation  and  the  flight  test  can  be  described  as  quite  promising,  the 
response  is  also  very  stable  with  the  exception  of  the  spiral  which  however  does  not 
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represent  a  serious  problem  since  it  has  a  time  to  double  of  about  20  seconds  However, 
the  pilot  should  be  aware  of  this  tendency  of  the  aircraft  to  avoid  trim  problems. 

B.  RECOMMENDATIONS 

It  is  recommended  that  the  stability  and  control  derivatives  be  calculated  for  other 
flight  conditions.  This  could  be  done  in  a  quick  and  accurate  way  using  the  same  programs 
in  Matlab.  Some  wind  tunnel  testing  would  also  verify  these  calculations  and  explore  some 
non-linear  flight  conditions  like  flight  at  high  angles  of  attack. 

The  purchase  of  a  portable  Unix  based  workstation  would  be  a  wise  step  towards 
the  integration  of  the  Matrixx/Systembuild  as  an  essential  tool  for  flight  testing  and  design. 
It  is  currently  necessary  to  transport  a  full  size  Sun  workstation  to  the  test  site.  Doing  so 
because  of  the  large  size  and  the  sensitivity  of  the  hardware  the  equipment  is  jeopardized 
during  the  transport  process. 

A  formal  nonlinear  parameter  estimation  approach  should  be  incorporated  into  the 
flight  test  program.  However  several  aspects  that  currently  are  not  known,  such  as  data 
acquisition  of  all  desired  information,  data  filtering  and  sensor  behavior  should  be 
addressed  prior  to  initiating  parameter  estimation. 

The  nonlinear  model  created  in  Systembuild  is  a  perfect  platform  for  designing 
control,  guidance  and  navigation  systems  which  will  be  implemented,  flight  tested  and 
improved  using  the  capabilities  of  Matrixx.  Stepping  through  the  rapid  prototype  design 
process  the  design  of  the  controller  is  significantly  simplified. 
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APPENDIX  A:  FILES  FOR  STABILITY  AND  CONTROL 

DERIVATIVES 


A.  FLIGHT  CONDITION 

%  File  for  Frog  data  which  change  with  flight  condition 

%  frogfcl.m 

%  Last  Update:  02  JAN  97 


g  =  32.174; 
Wlmg  =  30.42, 
Wrmg  =  32.13, 
Wng  =  5.18; 


%  Acceleration  due  to  gravity 
%  Weight  on  left  main  in  lbs 
%  Weight  on  right  main  in  lbs 
%  Weight  on  nose  gear  in  lbs 


Umph  =  60, 
Ufps  =  88.0; 
rho  =  .002327, 


%  Flight  speed  in  miles  per  hour 
%  Flight  speed  in  feet  per  second 
%  Air  density  in  slugs/(cubic  ft) 


Ixx=  12.52; 
Iyy  =  8.43, 
Izz=18  55, 
lxz  =  0: 


%  Moment  of  inertia  about  x-axis 
%  Moment  of  inertia  about  y-axis 
%  Moment  of  inertia  about  z-axis 
%  Assumed!!!!!!!!!!!!!!!!!!!!!!! 


LD=7, 


%  Lift  to  drag  ratio 


thetanaut  =0; 


%  Initial  pitch  angle 
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B.  FROG  DATA 

%  File  for  Frog  data  which  are  fixed 

%  frog.m 

%  Last  Update:  02  Mar  97 


ac  =  271; 

%  Aileron  chord  in  ft.# 

ai  =2.625, 

%  distance  from  centerline  to 

% 

inner  edge  of  aileron  in  ft.# 

alpha01  =  -2*pi/180, 

%  a.o.a.  for  zero  lift  (radians)# 

ao  =  4.417; 

%  distance  from  centerline  to 

% 

outer  edge  of  aileron  in  ft.# 

b=  10.58; 

%  Span  of  wing  in  ft# 

bt  =3.3125, 

%  Span  of  horizontal  tail  in  ft.# 

bv=1.25; 

%  Height  of  vertical  tail  in  ft.# 

cbar=  1.66, 

%  Mean  aerodynamic  chord  (m.a.c.) 

% 

in  feet# 

CLalphaafw  =5.87649; 

%  Lift  curve  slope  of  wing 

% 

airfoil  (NACA  2415)  in  per 

% 

radian# 

CLalphaaft  =5.72958; 

%  Lift  curve  slope  of  horizontal 

% 

tail  airfoil  (NACA  0010-34)  in  per 

% 

radian# 

CLalphaafV  =  2*pi; 

%  Lift  curve  slope  of  vertical 

% 

tail  airfoil  (flat  plate)  in  per 

% 

radian# 

CMac  =  -.045, 

%  Coefficient  of  moment  about 

% 

aero.  ctr.# 

ct  =0.968; 

%  m.a.chord  of  horizontal  tail  in  ft. 
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c4tail  =5.9086,  %5  8306, 

%  Location  of  quarter  chord  of 

% 

horizontal  tail  in  feet  from 

% 

proptip 

c4wing  =1.3108,  %1. 2483, 

%  Location  of  quarter  chord  of 

% 

wing  in  feet  from  proptip 

da0dde=625, 

%  Section  flap  effectiveness 

% 

for  33%  flap  (elevator) 

% 

Abbott  and  Doenhoff  p.  190 

da0ddr=675, 

%  Section  flap  effectiveness 

% 

for  38%  flap  (rudder) 

% 

Abbott  and  Doenhoff  p.  190 

deda  = .4, 

%  Downwash  angle  derivative 

% 

estimated  from  Perkins/Hage 

Df=  1.0416, 

%  Depth  of  fuselage  in  ft. 

e0  =  0, 

%  Assumed  epsilon  naught 

ee=8, 

%  Assumed  span  efficiency  factor 

g  =  32.174; 

%  gravitational  constant 

hac=241, 

%  Location  in  percent  chord  of 

% 

aero.  ctr.  (NACA  2415) 

it  =  (4.5+2)*pi/180, 

%  Incidence  angle  of  hor.  tail 

lewing  =  8958; 

%  Location  of  leading  edge  of  wing 

% 

in  feet  from  proptip 

letail  =5.667, 

%  Location  of  leading  edge  of 

% 

horizontal  tail  in  feet  from 

% 

proptip 

mg  =19.5/12, 

%  Location  of  main  gear  in  ft 

% 

from  proptip 

ng=-5/12, 

%  Location  of  nose  gear  in  ft 
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% 

from  proptip 

S=  17.5; 

%  Reference  (wing)  area  in  sq  ft 

Sr=.  501736; 

%  Rudder  area  in  sq.  ft 

St  =3.17448, 

%  Horizontal  tail  area  in  sq  ft  .# 

Sv  =0.98177, 

%  Vertical  tail  area  in  sq.  ft  # 

Wf=.75, 

%  Width  of  fuselage  in  ft.# 

ybar  =  b/4, 

%  Spanwise  location  of  mac. # 

zv=416, 

%  Vert,  tail  height  to  mac. 

% 

(estimated) 

Zwf  =  .5833; 

%  Vertical  height  of  wing 

% 

above  fuselage  C.L.  in  ft. 
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C.  PHYSICAL  CALCULATIONS 

%  File  to  calculate  physical  considerations 

%  physfrog.m 

%  Last  Update:  04  FEB  97 

%  Load  frog  data 

frog 

%  Load  flight  condition  data 

fcfrog 

%  Calculate  aircraft  weight 
W  =  Wlmg  +  Wrmg  +  Wng; 

%  Calculate  aircraft  mass 
m  =  W/g; 

%  Calculate  aspect  ratio  of  wing 
AR  =  bA2/S; 

%  Calculate  aspect  ratio  of  nor.  tail 
ARt  =  btA2/St, 

%  Calculate  aspect  ratio  of  vert,  tail 
ARv  =  bvA2/Sv, 

%  Calculate  longitudinal  center  of  gravity 

h  =  ((ng*Wng  +  mg*(Wlmg+Wrmg))/W-lewing)/cbar, 
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%  Calculate  "tail  length"  from  e.g.  to  horizontal  tail  a.c. 
%  same  for  horizontal  and  vertical 
It  =  c4tail  -  (lewing  +  h*cbar), 
lv  =  It; 

%  Calculate  "tail  length"  from  c/4  wing  to  c/4  tail 
ltprime  =  c4tail  -  c4wing, 

%  Calculate  hor.  tail  volume  coefficient 
VH  =  lt*St/(S*cbar); 

%  Calculate  vert,  tail  volume  coefficient  (yaw) 
W  =  lv*Sv/(b*S), 

%  Calculate  vert,  tail  volume  coefficient  (roll) 
Vprime  =  zv*Sv/(b*S); 

%  Unit  antisymmetrical  angle  of  attack  for  outer  and  inner 

%  edge  of  aileron  (See  Smetana  p.  141) 

antisymo  =  ao/(b/2);%  0.83497 

Cldatauo  =625; 

antisymi  =  ai/(b/2);%  0.49622 

Cldataui  =  .248, 

cacw  =  ac/cbar;%0. 16325 

tau  =  .48; 

%  for  yawing  moment  due  to  aileron,  see  p.  142,  Smetana 
eta  =  ai/(b/2),%0.49622  AR=6.3963 
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K  =  -.175; 

%  for  side  force  due  to  rudder  deflection,  see  Smetana  p.  145 
vratio  =  Sr/Sv;%  0.51 1052 
taur  =  .675; 

%  for  rolling  moment  due  to  sideslip,  See  Raymer,  Fig.  16.21,  p.  439 
ClBwCL  =  -.04, 
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D.  NONDEVfENSIONAL  DERIVATIVES  CALCULATIONS 

%  File  to  calculate  nondimensional  derivatives 

%  ndderiv.m 

%  Last  Update:    12  FEB  97 

%  Load  Frog  data  with  flight  condition 
phfrog 

%  Calculate  coefficients  of  lift  and  drag 
CL  =  W/(.5*rho*UfpsA2*S), 
CD  =  CL/LD, 
D=CD*(.5*rho*UfpsA2*S), 

%  Calculate  lift  curve  slope  of  wing  in  per  radian 
CLalphaw  =  CLalphaafw/(l+CLalphaafw/(pi*ee*AR)), 

%  Calculate  lift  curve  slope  of  horizontal  tail  in  per  radian 
CLalphat  =  CLalphaaft/(l+CLalphaan7(pi*ee*ARt)); 

%  Calculate  lift  curve  slope  of  vertical  tail  in  per  radian 
CLalphav  =  CLalphaafv/(l+CLalphaafV/(pi*ee*ARv)), 

%  Calculate  change  in  hor.  tail  lift  with  change  in  elevator 
dcLtdde  =  daOdde  *  CLalphat,  %  per  radian 

%  Calculate  change  in  vert,  tail  lift  with  change  in  rudder 
dcLvddr  =  daOddr  *  CLalphav,  %  per  radian 
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%  Calculate  zero  lift  pitching  moment 
CMO  =  CMac  +  VH  *  CLalphat  *  (it  +  eO), 

%  Calculate  CMalpha  in  per  radian 

CMalpha  =  CLalphaw*((h-hac)-VH*(CLalphat/CLalphaw)*(l-deda)), 

%  Calculate  change  in  aircraft  lift  with  change  in  elevator 
CLdelta  =  dcLtdde*(St/S),      %  per  radian 

%  Calculate  change  in  aircraft  pitching  moment  w.  chng  in  elevator 
CMde  =  - 1  *  VH*dcLtdde;       %  per  radian 

%  Calculate  angle  of  attack  and  elevator  angle  for  trimmed  flight 

% 

%         CM  =  CMO  +  CMalpha*alpha  +  CMde*de 

%         CI  =  CLalphaw*  alpha  +  CLdelta*de 

% 

% 


%  |  CLalphaw    CLdelta  |  |  alpha  |      |  CL    | 

%  |  111  =  11 

%  |_CMalpha      CMde |  |_de_|      |_-CM0 

% 

%  A  *      X       =      C 

°/o 


0 

A  =  [  CLalphaw  CLdelta 
CMalpha     CMde    }; 
C  =  [     CL 

-1*CM0]; 


57 


X  =  inv(A)*C, 

atrim  =  X(l,l)  %  trim  a.o.a.  in  radians 

etrim  =  X(2, 1 )  %  trim  elevator  in  radians 

%  Calculate  change  in  yawing  moment  with  change  in  rudder 

%  "rudder  power" 

%  assumes  VF/Vinfinity  =  1 

Cndr  =  - 1  *  W*dcLvddr,         %  in  per  radian 

%  Calculate  CnB  contribution  from  vert,  tail 

%  CnB  =  CLalphav*W*(W/Vinfii%)A2*(l-dsigma/dbeta) 

%  assumes  W/Vinfinity  =  1  and  dsigma/dbeta  =  0 

CnB  =  CLalphav*  W;%  in  per  radian 

%  Calculate  change  in  rolling  moment  with  change  in  sideslip 

%  First  calculate  dihedral  contribution  from  wing 
%  Raymer  p.  439 

ClBwf  =  -1.2*sqrt(AR)*Zwf*(Df+Wf)/bA2, 
ClBw  =  ClBwCL*CL+ClBwf, 

%  Next  calculate  contribution  from  fin 

%  ClBv  =  -l*Clalphav*Vprime*(VFA^infinity)A2*(l-dsigma/dbeta) 

%  Assume  W/Vinfinity  =  1  and  dsigma/dbeta  =  0 

ClBv  =  - 1  *CLalphav* Vprime,  %  in  per  rad 

%  Combine  ClBg  and  ClBv  into  C1B 
C1B  =  ClBw  +  ClBv,   %  in  per  rad 
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%  Calculate  "aileron  power",  Clda 
%  See  Smetana  pp.  139-141 
Cldatau  =  Cldatauo  -  Cldataui, 
Clda  =  Cldatau*tau,     %  in  per  radian 

%  Calculate  change  in  yawing  moment  w.  aileron  deflection 
Cnda  =  2*K*CL*Clda,  %  in  per  radian 

%  Calculate  side  force  due  to  yaw 

%  By  Smetana  p.  107 

CyB  =  -.31,      %  in  per  radian 

%  Calculate  side  force  due  to  rudder 

Cydr  =  CLalphav*taur*Sv/S,  %  in  per  radian 

%  Calculate  side  force  due  to  aileron 
%  By  Smetana,  p.  138 
Cyda  =  0, 

%  Calculate  rolling  moment  due  to  rudder 
Cldr  =  Cydr*zv/b,       %  in  per  radian 

%  Calculate  change  in  drag  due  to  change  in  elevator 

%  Smetana  pp  95-100 

%  Using  Figure  25  at  6  degrees  aoa 

CDde  =  ((.  16-.03)/(20*pi/180))*St/S;  %  in  per  radian 
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%  Calculate  change  in  drag  with  change  in  aoa 

%  Smetana  pp.  64-65 

%  Assuming  dCDO/dalpha  is  negligible 

CDalpha  =  2*CL*CLalphaw/(pi*ee*AR),      %  in  per  radian 

%  Calculate  change  in  pitching  moment  w.r.t  alphadot 
%  Smetana  pp.  78-81,  etat  assumed  =  1 
CMalphadot  =  -2*CLalphat*deda*(ltprime/cbar)* ... 
(lt/cbar)*(St/S);  %  in  per  radian 

%  Calculate  change  in  lift  with  pitch  rate 

%  Smetana  pp.  82-85 

%  Neglecting  wing  contribution,  assuming  etat  =  1 

CLq  =  2*(lt/cbar)*CLalphat*(St/S);  %  in  per  radian 

%  Calculate  change  in  lift  with  alphadot 

%  Smetana  pp.  75-76 

CLalphadot  =  -1  *CMalphadot/(lt/cbar);         %  in  per  radian 

%  Calculate  change  in  pitching  moment  w.  pitch  rate 

%  Smetana  pp.  87-88 

%  Assuming  etat  =  1 

CMq  =  -2*(cbar/4-h*cbar)*abs(cbar/4-h*cbar)*CLalphaw/(cbarA2) 

2*(lt/cbar)A2*CLalphat*(St/S),  %  in  per  radian 

%  Calculate  roll  damping 

%  Smetana  pp.  122-125  %  Clp(a0:2pi)=-0.475 

%  Neglecting  contribution  from  vertical  tail 
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Clp  =  -.475*(AR+4)/(2*pi*AR/CLalphaw+4);  %  in  per  radian 

%  Calculate  change  in  yawing  moment  due  to  rolling 

%  Smetanapp  126-129 

%  Neglecting  contribution  from  vertical  tail 

Cnp  =  -1  *CL/8,  %  in  per  radian 

%  Calculate  change  in  side  force  with  yaw  rate 

%  From  Schmidt  p.  3-23 

%  Assume  etat  =  1 

Cyr  =  2*  W*CLalphav,  %  in  per  radian 

%  Calculate  change  in  rolling  moment  w.  yaw  rate 

%  Schmidt  p.  3-24 

%  Tail  contribution 

Clrv  =  (zv/b)*Cyr;       %  in  per  radian 

%  Wing  contribution 

Clrw  =  CL/4;   %  in  per  radian 

%  Total 

Clr  =  Clrv  +  Clrw,       %  in  per  radian 

%  Calculate  yaw  damping 

%  Schmidt  p.  3-25 

%  Tail  contribution 

Cnrv  =  -1  *(rv/b)*Cyr,  %  in  per  radian 

%  Wing  contribution  from  Smetana  p.  136 

CDO  =  CD-CLA2/(pi*ee*AR); 

Cnrw  =  -.02*CLA2-.3*CD0,  %  in  per  radian 
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%  Total 

Cnr  =  Cnrv  +  Cnrw;    %  in  per  radian 

%  The  following  3  derivatives  are  negligible  and  taken  to  be  0 
CDq  =  0,  %  in  per  radian 

Cyq  =  0;  %  in  per  radian 

Cyp  =  0,  %  in  per  radian 

%  A  few  misc.  calculations 

%  Static  Margin/Neutral  Point 
statmar  =  CMalpha/(-l  *CLalphaw), 
hn  =  statmar  +  h 
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E.  DIMENSIONAL  DERIVATIVES  CALCULATIONS 

%  File  to  calculate  dimensional  derivatives 

%  frogdder.m 

%  Last  Update:    12  FEB  97 

%  Run  nondimensional  derivative  program 
ndfrog 

%  Calculate  dynamic  pressure 

qbar  -  5*rho*UfpsA2;  %  ft  lbs 

Malpha  =  CMalpha*qbar*S*cbar/Iyy,  %  per  secondA2 

Mq  =  CMq*(cbar/(2*Ufps))*qbar*S*cbar/Iyy, 

%  per  second 

Malphadot  =  CMalphadot*(cbar/(2*Ufps))*qbar*S*cbar/Iyy, 

%  per  second 

Xu  =  -2*CD*qbar*S/(m*Ufps),  %  per  second 

Zu  =  -2*CL*qbar*S/(m*Ufps);  %  per  second 

Zalphadot  =  CLalphadot*(cbar/(2*Ufps))*(qbar*S/m), 

%  ft  per  second 

Zq  =  CLq*(cbar/(2*Ufps))*(qbar*S/m),        %  ft  per  second 
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Mu  =  0;  %  per  ft  second 

Xde  =  -1  *CDde*qbar*S/m;  %  ft  per  secondA2 

Zde  =  - 1  *CLdelta*qbar* S/m,  %  ft  per  secondA2 

Mde  =  CMde*qbar*S*cbar/Iyy;  %  per  second A2 

Xalpha  =  (CL  -  CDalpha)*qbar*S/m,  %  ft  per  secondA2 
Zalpha  =  - 1  * (CLalphaw+CD) * qbar*  S/m;       %  ft  per  secondA2 

YB  =  CyB*qbar* S/m,  %  ft/secondA2 

LB  =  ClB*qbar*S*b/Ixx;  %  l/secondA2 

NB  =  CnB  *qbar*  S  *b/Izz,  %  1  /secondA2 

Yp  =  Cyp*b*qbar*S/(2*Ufps*m);  %  ft/sec 

Yr  =  Cyr*b*qbar*S/(2*Ufps*m);  %  ft/sec 

Lp  =  Clp*(b/(2*Ufps))*qbar*S*b/Ixx,  %  1/sec 

Np  =  Cnp*(b/(2*Ufps))*qbar*S*b/Izz;  %  1/sec 
Lr  =  Clr*(b/(2*Ufps))*qbar*S*b/Ixx;%  1/sec 
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Nr  =  Cnr*(b/(2*Ufps))*qbar*S*b/Izz,  %  1/sec 

Ydr  =  -1  *Cydr*qbar*S/m,  %  ft/secA2 

Yda  =  0,  %  ft/secA2 

Ldr  =  Cldr*qbar*S*b/Ixx,  %  l/secA2 

Lda  =  Clda*qbar*S*b/Ixx;  %  1/sec  A2 

Ndr  =  Cndr*qbar*S*b/Izz;  %  l/secA2 

Nda  =  Cnda*qbar*S*b/Izz,  %  l/secA2 

Malphaprime  =  Malpha  +  Malphadot*(Zalpha/Ufps), 

Mqprime  =  Mq  +  Malphadot; 

Mdeprime  =  Mde  +  Malphadot  *(Zde/Ufps), 
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F.   STABILITY  AND  CONTROL  DERIVATIVES 

%  File  to  get  values  of  dimensional  and  nondimensional  derivatives 
%  Last  Update.  04  FEB  97 

%  Run  programs  to  calculate  derivatives 
ddfrog 

%  Nondimensional  Derivatives 

CL 

CD 

CDO 

CLalphaw 

CMalpha 

CDalpha 

CLalphadot 

CMalphadot 

CLq 

CMq 

CLdelta 

CDde 

CMde 

CyB 

CnB 

C1B 

Cyp 

Cnp 

Clp 

Cyr 
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Cnr 

Clr 

Cyda 

Cnda 

Clda 

Cydr 

Cndr 

Cldr 

CDq 

Cyq 

%  Longitudinal  Dimensional  Derivatives 

Xu 

Xalpha 

Xde 

Zu 

Zalpha 

Zalphadot 

Zq 

Zde 

Mu 

Malpha 

Malphadot 

Mq 

Mde 

%  Lateral/Directional  Dimensional  Derivatives 
YB 
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Yp 

Yr 

Yda 

Ydr 

LB 

Lp 

Lr 

Lda 

Ldr 

NB 

Np 

Nr 

Nda 

Ndr 
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APPENDIX  B:  MOMENTS  OF  INERTIA  CALCULATIONS 


For  the  formula: 


_  Wm+1<»(,Lshort+Lu>ng)¥m+s  Jl  wmzm        2coxr3  ,r3         x 

FM+S  g  J^K^  SHORT   +l-LONG) 


4nl 


the  following  data  were  used  for  the  calculation  of  the 

Constant  values: 

Weight  of  the  model, 
Weight  of  the  unit  length  chain, 
Length  of  short  chain, 
Length  of  long  chain, 
Gravitational  constant 


moments  of  inertia 


WM  =  67.73  lbs. 
(o  =  0.061887  lbs/ft. 
L  SHORT  =  13.5//. 
Lwng  =  15//. 
g  =  32. 1472  ftls2. 

(adjusted  for  latitude  and  elevation) 


Variable  values: 
I xx 


<YY- 


fzz 


Distance  from  pivot  to  center  of  gravity 

of  model  and  support, 

Distance  from  pivot  to  center  of  gravity 

of  model, 

Average  period  of  model  and  support, 
Distance  from  pivot  to  center  of  gravity 
of  model  and  support, 

Distance  from  pivot  to  center  of  gravity 

of  model, 

Average  period  of  model  and  support, 
Distance  from  pivot  to  center  of  gravity 
of  model  and  support, 

Distance  from  pivot  to  center  of  gravity 

of  model, 

Average  period  of  model  and  support, 


ZM+s=  13.91166// 
ZM  =  14  ft. 

Pm+s  =  4. 1847 ^  sec 
ZM+s  =13.91166  ft. 

ZM=\4ft. 

/W  =  4. 16475  sec 

^+5=13.41775^. 

Za/=13.5//. 
/V+s  =  4.1455  sec. 
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APPENDIX  C:  MATLAB  FILES  OF  EQUATIONS  OF  MOTION 


A.  ABMAT.M 

%  File  for  modelling  the  nonlinear  equations  of  motion 

%  Last  Update:  21  JAN  1997 

function  accel  =  abmat(x) 

%  calculates  the  accelerations  (angular  and  linear)  due  to 

%  aero  forces,  w  X  v;  gravity. 

%  Variables  brought  from  workspace: 

%  x  =  combination  vector  [contrl  inputs,  state  variables,  euler  angles] 

%  (da,  de,  dr,  dtrt,  u,  v,  w,  p,  q,  r,  phi,  theta,  psi) 

%  Variables  called  from  function  "frogdata" 

%  rho  =  air  density 

%  b  =  wing  span 

%  c  =  wing  cord 

%  s  =  wing  area 

%  Cfo  =  Steady  state  force  term 

%  Cfu  =  Stability  derivative  for  control  inputs 

%  m  =  airplane  mass 

%  lb  =  inertia  tensor  matrix  (body  frame) 

%  To  =  Thrust  scale  term 

%  Pe  =  Engine  position  matrix 

%  get  the  aircraft  data 
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[uO,wO,rho,Cfx,Cfo,Cm,Cfxdot,s,b,c,m,Pe  Jo,Ib]  =  frog_data, 
%  seperate  the  combined  vector  into  seperate  elements 

u  =  [x(l);x(2);x(3)]; 

dm  =  x(4), 

state  =  [x(5),  x(6);  x(7);  x(8);  x(9);  x(10)], 

lambda  =  [x(l  1),  x(12);  x(13)]; 

%%%%%%  calculate  velocity  wrt  the  airmass  and  form  state  vector 
%%%%%%  that  will  be  used  to  calculate  the  aerodynamic  forces/moments 

statel  =  [state(l)-uO;  state(2);  state(3)-w0;  x(8);  x(9);  x(10)]; 

%%%%%%  calculate  total  velocity,  vt 

vt  =  (state(l)A2  +  state(2)A2  +  state(3)A2)A.5; 

%  calculate  qbar 

qbar=.5*rho*(vtA2), 

%  calculate  Ml 

Ml  =  diag([l/vt,  1/vt,  1/vt,  (.5*b)/vt,  (.5*c)/vt,  (5*b)/vt]); 

%  calculate  M2 

M2  =  diag([0,  0,  (,5*c)/(vtA2),  0,  0,  0]), 
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%  calculate  Sprime 

Sprime  =  diag([-l,  1,  -1,  b,  c,  b]*s); 

%  calculate  Mu 

Mu  =  inv([eye(3)*m,zeros(3);zeros(3),Ib]); 

%  calculate  Tw2b 

alpha  =  state(3)/vt, 

beta  =  state(2)/vt; 

Tl  =  [cos(alpha),  0,  sin(alpha),  0,1,0,  -sin(alpha),  0,  cos(alpha)], 

T2  =  [cos(beta),  -sin(beta),  0,  sin(beta),  cos(beta),  0,  0,0,1]; 

Tw2b  =  [Tr*T2',  zeros(3),  zeros(3),  T1'*T2']; 

%  calculate  Chi 

Chi  =  eye(6)  -  Mu*Tw2b*qbar*Sprime*Cfxdot*M2, 

%  calculate  Propulsion  matrix 

Prop  =  [  eye(3); 
0,-Pe(3),Pe(2), 
Pe(3),0,-Pe(l), 
-Pe(2),Pe(l),0]; 
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%  calculate  gravity  vector  and  rotation  matrix 

Rot  =  [1,  0,  -sin(lambda(2)), 

0,  cos(lambda(l)),  cos(lambda(2))*sin(lambda(l)); 

0,  -sin(lambda(l)),  cos(lambda(2))*cos(lambda(l))]; 

Ru2b  =  [Rot,zeros(3)], 

g=[0;0;  32.174]; 
FgU  =  m*g, 

%  put  the  components  due  to  gravity,  thrust,  and  control  surface  deflections 
%  together  for  their  contribution  to  x-dot 

thrust  =  Prop*To*dtrt; 
gravity  =  Ru2b*FgU; 
ctrl=qbar*(Tw2b*(Sprime*(Cfo  +  (Cfu*u)))); 

xdotu=(Mu  *  (ctrl+thrust+gravity)) ; 

%  calculate  rotation  matrix 

omegax  =  [0,-state(6),state(5),state(6),0,-state(4),-state(5),state(4),0], 
wxlb  =  (-inv(Ib))*(omegax*Ib), 

Rot  =  [-omegax,  zeros(3),  zeros(3),  wxlb], 

%  rotation  component  of  xdot  (w  X  v) 
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xdotrot  =  Rot*  state, 

%  state  vector  feedback  component  xdot 

xdotcfx  =qbar*(Mu*(Tw2b*(Sprime*(Cfx*(Ml  *statel))))); 

%  add  three  components  of  xdot  together  and  premult  by  inv(Chi) 

xdot=  inv(Chi)  *  (xdotrot+xdotcfx+xdotu) , 

%  calc  accel  that  a  strap-down  IMU  would  measure 

xdotb=xdot-xdotrot-Ru2b*g, 

%  put  together  for  the  return  vector 

%accel=[xdotb(  1  );xdotb(2),xdotb(3),xdot], 

%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%%0/o%%%0/o0/o%%%%%%% 

%%%%%% 

%  return  xdot  only 

accel=xdot, 
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B.  FROG_DAT.M 

%  File  for  inputting  all  aircraft  data  for  equations  of  motion 

%  Last  Update:  21  JAN  1997 

function  [uO,wO,rho,Cfx,Cfo,Cfu,Cfxdot,s,b,c,m,Pe,To,Ib]  =frog_data 

uO  =  88, 
w0  =  0, 

%  Sea  level-  std  day 

rho  =  .0023769, 

%  derivative  matrix  due  to  state  variables 
%  rows:  [CD  CY  CL  CI  Cm  Cn] 

%  col:  [u  v  w  p  q  r] 

Cfx  =  [0    0.23  0  0  0, 
0-31000.1151, 
0    0  4.3034  0  3.35  0, 
0 -.0509  0-.3702  0.il  19; 

0  0-0.5565  0-8.8818  0; 
0  .0575  0  -.0537  0  -.0669]; 

%  derivative  matrix  due  to  control  inputs 
%  rows:  [same  as  above] 
%  col:  [elev  rud  ail] 
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Cfu  =  [0676  0  0, 
0  0926  0; 
.3914  0  0; 
0.0036.1810, 
-1.0469  0  0, 
0  -.0388  -.0272], 

%  derivative  matrix  due  to  x-dot 

Cfxdot  =  [000000, 

0  0  0  0  0  0, 

0  0  1.3877  0  0  0, 

0  0  0  0  0  0, 

0  0-3.7115000, 

0  0  0  0  0  0]; 

%  steady  state  force  vector 

Cfo  =  [0614,  0,  .4295,  0;  0,  0], 

%  physical  dim. 
%WT  =67.73  LBS. 

m  =  2.1051; 
s=17.5, 
b=  10.58, 
c=  1.66; 
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%  engine  data 

Pe  =  [17.63/12;  0,-14.92/12], 
To  =  [9.6757  ,0,0], 

%  inertia  tensor  matrix 

lb  =  [  12.52  0  0,  0  8.43  0,  0  0  18.55]; 
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APPENDIX  D:  NUMERICAL  LINEARIZATION  RESULTS 


%  Numerical  results  from  linearization  of  equations  of  motion 

%  Last  Update:  20  FEB  97 

%  The  state  vector  at  trim  will  be: 

x=[  87.9999;  0.0000,  0.1593;  0.0000;  0.0000,  0.0000,  0  0000;  0.0018,  0.0000]' 
%  the  inputs  at  trim  are: 

u=[ -0.0431;  0.0000,  0.0000,  0.9805]' 
%  while  the  outputs  are  as  defined: 

y  =[0.0000,  88.0000]' 

%  Then  the  matrices  of  the  linear  model  are  the  following: 

A=[  -0.1014  0.0000  0.1722  0    -0.1532          0     0.0000      -32.1739  0 

0.0000  -0.2183  0.0000  0.1593    0      -87.4705  32.1739    0.0000  0 

-0.7162  0.0000  -3.7510  0     84.6195         0     -0.0002    -0.0578  0 

0.0000  -0.0682  0.0000  -3.0280  0.0000     0.9165  0.0000    0.0000  0 

0.0412  0.0000  -0.1532  0      -3.7244         0       0.0000     0.0007  0 

0.0000  0.0599  0.0000  -0.3002  0.0000  -0.3683  0.0000     0.0000  0 

0  0  0  1.0000  0.0000     0.0018  0.0000     0.0000  0 

0  0  0  0          1.0000     0.0000     0.0000      0  0 

0  0  0  0          0.0000     1.0000     0.0000    0.0000  0] 
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B  =[-5.1184    0    0  4.5963 

0.0000  7.0847   0    0 
-29.6178    0    0    0 

0  0000  0.4995  24.6412  0 
-32.8288    0    0  -1.4271 

0.0000  -3.5636  -2.4685  0 
0      0    0    0 
0      0    0    0 
0      0    0    0] 


C=[  0.0000    0  -0.0114    0    0    0 
1.0000  0.0000  0.0018    0    0    0 


0       1.000       0 
0       0  0] 


D=[0      0      0      0,0      0      0      0] 


%  The  eigenvalues  of  the  A  matrix  are: 


Eigenvalue            Damping 

Freq.  (rad/sec) 

0.0474 

1.0000 

0.0474 

0 

■1.0000 

0 

-0.0293  +  0.5597i 

0.0522 

0.5604 

-0.0293  -  0.5597i 

0.0522 

0.5604 

-0.1964 +  2.4972i 

0.0784 

2.5049 

-0.1964-2.4972i 

0.0784 

2.5049 

-3.2691 

1.0000 

3.2691 

-3.7591  +3.5964i 

0.7226 

5.2024 

-3.7591- 3. 5964i 

0.7226 

5.2024 
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E.  ADDITIONAL  SUPERBLOCK  DIAGRAMS 
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Fig.  E.l  Frog  superblock  diagram 
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Fig.  E.2  Integrators  Superblock 
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APPENDIX  F:  FLIGHT  TEST  RESULTS 


60 

40 

20 

0 

Asimulation         j 

flight  test]     \l 

K                 1 

-40 
-60 

15 


Time-seconds 


Fig.  F.l  Short  Period  Response-  Run  #1 
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Fig.  F.2  Short  Period  Response-  Run  #2 
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Fig.  F.3  Short  Period  Response-  Run  #3 
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Fig.  F.4  Short  Period  Response-  Run  #4 
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Fig.  F.5  Phugoid  Response-  Run  #1 
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Fig.  F.6  Phugoid  Response-  Run  #2 
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Fig.  F.7  Phugoid  Response-  Run  #3 
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Fig.  F.8  Steady  Sideslip-  Run  #1 
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Fig.  F.9  Steady  Sideslip-  Run  #2 
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Fig.  F.10  Steady  Sideslip-  Run  #3 
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Fig.  F.ll  Steady  Sideslip-  Run  #  4 
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Fig.  F.12  Steady  Sideslip-  Run  #5 
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Fig.  F.13  Dutch  Roll/Spiral  Response-  Run  #1 
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Fig.  F.14  Dutch  Roll/Spiral  Response-  Run  #2 
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Fig.  F.15  Dutch  Roll/Spiral  Response-  Run  #3 
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Fig.  F.16  Dutch  Roll/Spiral  Response-  Run  #4 
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